//Includes
#include <stdio.h>
#include <unistd.h>
#include "stdlib.h"
#include "math.h"
#include "string.h"
#include "sys/time.h"
#include <mpi.h>

//Defines
#define NUM_ARGS  7
#define WAIT_TIME 5
#define SOURCE        	   0

//Prototypes
double* allocateArray(const int reducedRow, const int reducedColumn);
void initArray(double *array, const int numRow, const int numColumn, const int reducedRow, const int reducedColumn);
void freeArray(double *array);

double** allocateMatrix(const int numRows, const int numColumns);
double** convertArrayToMatrix(double *array, const int numRow, const int numColumn,
		 	 	 	 	 	 const int reducedRow, const int reducedColumn,
		 	 	 	 	 	 const int rowIsOdd, const int columnIsOdd);
void freeMatrix(double **matrix, const int numRows);

void printMatrix(double **matrix, const int numRows, const int numColumns, int isReversed);

void calculateEulerFormula(double *array, const int reducedRow, const int reducedColumn, const int numTimeStep,
						   const double discredTime, const double subdivisionSize,
						   const int rowIsOdd, const int columnIsOdd);

void calculateEulerFormulaParallel(double *array, const int reducedRow, const int reducedColumn, const int numTimeStep,
						   	   	   const double discredTime, const double subdivisionSize,
						   	   	   const int rowIsOdd, const int columnIsOdd,
						   	   	   int *aggregationArray, int rank, int numtasks);

int* aggregate(int numtasks, int reducedRow, int reducedColumn);
void updateNeighbor(double *previousArray,int *aggregationArray, int rank, int numtasks, int reducedColumn);

int main(int argc, char *argv[])
{
	if(argc == NUM_ARGS)
	{
		MPI_Init(&argc,&argv);
		int numtasks, rank;
		MPI_Comm_size(MPI_COMM_WORLD, &numtasks);
		MPI_Comm_rank(MPI_COMM_WORLD, &rank);


		double *array = NULL;
		double **matrix = NULL;
		int reducedRow, reducedColumn, nbReducedElement;
		int rowIsOdd, columnIsOdd;
		int isReversed = 0;


		int numRows = atoi(argv[1]);           // n
		int numColumns = atoi(argv[2]);        // m
		int numTimeStep = atoi(argv[3]);       // np
		double discredTime = atof(argv[4]);     // td
		double subdivisionSize = atof(argv[5]); // h
		int nbProcessors = atoi(argv[6]);      // nbproc

		// minuteur
		double timeStart, timeEnd, TexecParallel, TexecSequential;
		struct timeval tp;

		// parameter validation
		if((numRows < 3) || (numColumns < 3) || (numTimeStep < 1) || (discredTime <= 0.0f) ||
		   (subdivisionSize <= 0.0f) || (nbProcessors < 1))
		{
			if(rank == SOURCE)
				printf("Arguments are invalid.\n");
			MPI_Finalize();
			exit(1);
		}

		// we make sure the matrix is vertical...if not we rotate it
		if (numColumns > numRows)
		{
			numRows = atoi(argv[2]);
			numColumns = atoi(argv[1]);
			isReversed = 1;
			if(rank == SOURCE)
				printf("m > n :: the matrix is horizontal, it will be reversed and treated as a vertical\n");
		}
		else if (rank == SOURCE)
			printf("m < n :: the matrix is vertical\n");


		// measures of the reduced matrix (perimeter removed and we only keep top left quadrant - we will use symetry to complete)
		reducedRow = (int) (((numRows - 1) / 2));
		reducedColumn = (int) (((numColumns - 1) / 2));
		nbReducedElement = reducedRow * reducedColumn;
		rowIsOdd = (numRows % 2);
		columnIsOdd = (numColumns % 2);


		// validate that the number of row to resolve is even or higher than the number of processor used.
		// those are the minimum required condition for the parallel algorithm (aggregation/communication) to work
		if (numtasks > reducedRow)
		{
			if(rank == SOURCE)
				printf("the amount of processor used %d is too high for such a small plate. Please reduced to %d or increase size of the plate and start application again", numtasks, reducedRow);
			MPI_Finalize();
			exit(1);
		}


		// the reduced matrix will be stored in an array
		array = allocateArray(reducedRow, reducedColumn);

		// ************************** Parallel version start **************************
		int *aggregationArray = NULL;

		initArray(array, numRows, numColumns, reducedRow, reducedColumn); //initialise value of the matrix
		aggregationArray = aggregate(numtasks, reducedRow, reducedColumn);//specify the aggregation

		if(rank == SOURCE)
		{
			printf("************************** Parallel version **************************\n");
			printf("numtask: %d \n", numtasks);
			printf("reducedRow: %d \n", reducedRow);
			printf("reducedColumn: %d \n\n", reducedColumn);

			printf("processor: index of 1st element to treat \n");
			int i;
			for(i = 0; i <numtasks; i++)
				printf("%d: %d \n", i, aggregationArray[i]);

			// Print initialised matrix
			matrix = convertArrayToMatrix(array, numRows, numColumns,  reducedRow, reducedColumn, rowIsOdd, columnIsOdd);
			printf("Matrix after initialisation:\n");
			printMatrix(matrix, numRows, numColumns, isReversed);

			// Start Timer
			gettimeofday (&tp, NULL);
			timeStart = (double) (tp.tv_sec) + (double) (tp.tv_usec) / 1e6;
		}

		// all the work is done here (computing/communication/data gathering)
		calculateEulerFormulaParallel(array, reducedRow, reducedColumn, numTimeStep, discredTime, subdivisionSize, rowIsOdd, columnIsOdd, aggregationArray, rank, numtasks);

		if(rank == SOURCE)
		{
			// End timer
			gettimeofday (&tp, NULL);
			timeEnd = (double) (tp.tv_sec) + (double) (tp.tv_usec) / 1e6;
			TexecParallel = timeEnd - timeStart; //Temps d'execution en secondes
			printf("Execution time: %f sec.\n", TexecParallel);

			// Print Matrix
			printf("Matrix after Euler:\n");
			matrix = convertArrayToMatrix(array, numRows, numColumns,  reducedRow, reducedColumn, rowIsOdd, columnIsOdd);
			printMatrix(matrix, numRows, numColumns, isReversed);
			freeMatrix(matrix, numRows);
		}

		if (rank != SOURCE)
			freeArray(array);

		// ************************** End of Parallel Version **************************


		// ************************** Sequential version start **************************
		if(rank == SOURCE)
			{
			printf("************************** Sequential version start **************************\n");

			// initialize matrix
			initArray(array, numRows, numColumns, reducedRow, reducedColumn);//initialise value of the matrix

			// Print matrix
			matrix = convertArrayToMatrix(array, numRows, numColumns,  reducedRow, reducedColumn, rowIsOdd, columnIsOdd);
			printf("Matrix after initialisation:\n");
			printMatrix(matrix, numRows, numColumns, isReversed);

			// Start Timer
			gettimeofday (&tp, NULL);
			timeStart = (double) (tp.tv_sec) + (double) (tp.tv_usec) / 1e6;

			// Compute
			calculateEulerFormula(array, reducedRow, reducedColumn, numTimeStep, discredTime, subdivisionSize, rowIsOdd, columnIsOdd);

			// End timer
			gettimeofday (&tp, NULL);
			timeEnd = (double) (tp.tv_sec) + (double) (tp.tv_usec) / 1e6;
			TexecSequential = timeEnd - timeStart; //Temps d'execution en secondes
			printf("Execution time: %f sec.\n", TexecSequential);

			// Print Matrix
			printf("Matrix after Euler:\n");
			matrix = convertArrayToMatrix(array, numRows, numColumns,  reducedRow, reducedColumn, rowIsOdd, columnIsOdd);
			printMatrix(matrix, numRows, numColumns, isReversed);

			freeMatrix(matrix, numRows);
			freeArray(array);

			// print execution times, number of processor used and acceleration achieved
			printf("Parallel execution time: %f sec.\n", TexecParallel);
			printf("Sequential execution time: %f sec.\n", TexecSequential);
			printf("Acceleration: %#1.2f \n", TexecSequential/TexecParallel);
			printf("Number of processor used: %d \n", numtasks);
			}
		// ************************** End of Sequential Version **************************

		MPI_Finalize();
	}
	else
	{
		printf("The number of arguments is invalid.\n");
	}

	return 0;
}

void updateNeighbor(double *previousArray, int *aggregationArray, int rank, int numtasks, int reducedColumn)
{
	MPI_Status statusFromLeft;
	MPI_Status statusFromRight;

	// send to left neighbor
	if(rank > SOURCE)
	{
		MPI_Send(&previousArray[aggregationArray[rank]], reducedColumn, MPI_DOUBLE, rank - 1, 0, MPI_COMM_WORLD);
	}

	//Receive message from Right neighbor
	if(rank < numtasks-1)
	{
		MPI_Recv(&previousArray[aggregationArray[rank+1]], reducedColumn, MPI_DOUBLE, rank + 1, 0, MPI_COMM_WORLD, &statusFromRight);
	}

	// send to right neighbor
	if(rank < numtasks-1)
	{
		MPI_Send(&previousArray[aggregationArray[rank+1]-reducedColumn], reducedColumn, MPI_DOUBLE, rank + 1, 0, MPI_COMM_WORLD);
	}

	//Receive message from Left neighbor
	if(rank > SOURCE)
	{
		MPI_Recv(&previousArray[aggregationArray[rank]-reducedColumn], reducedColumn, MPI_DOUBLE, rank - 1, 0, MPI_COMM_WORLD, &statusFromLeft);
	}
}

// this function return an array with the 1st element of the reduced matrix that each processor(index) must treat
// every processor treat an even number of element except the last one who may treat less
int* aggregate(int numtasks, int reducedRow, int reducedColumn)
{
	int *array;
	int i;
	int nbElement = reducedColumn * reducedRow;
	int avgElementPerProc = ceil((float)nbElement / numtasks);

	array = (int*) malloc((numtasks) * sizeof(int));

    if(array == NULL)
    {
        printf("Out of memory\n");
        exit(1);
    }

    for(i = 0; i < numtasks; i++)
    {
    	array[i] = i * avgElementPerProc;
    }

	return array;
}

double* allocateArray(const int reducedRow, const int reducedColumn)
{
	double *array;

    array = (double*) malloc((reducedRow * reducedColumn) * sizeof(double));

    if(array == NULL)
    {
    	printf("Out of memory\n");
    	exit(1);
    }

    return array;
}

double** allocateMatrix(const int numRows, const int numColumns)
{
	double **matrix;
    int rowCounter;

    matrix = (double**) malloc(numRows * sizeof(double*));

    if(matrix == NULL)
    {
        printf("Out of memory\n");
        exit(1);
    }

    for(rowCounter = 0; rowCounter < numRows; rowCounter++)
    {
        matrix[rowCounter] = (double*) malloc(numColumns * sizeof(double));

        if(matrix[rowCounter] == NULL)
        {
                printf("Out of memory\n");
                exit(1);
        }
    }

    return matrix;
}

void freeMatrix(double **matrix, const int numRows)
{
        int rowCounter;

        for(rowCounter = 0; rowCounter < numRows; rowCounter++)
        {
                free( matrix[rowCounter] );
        }

        free( matrix );
}

// Function that return a the full matrix that we rebuild with data contained in a reduced one
double** convertArrayToMatrix(double *array, const int numRow, const int numColumn,
							 const int reducedRow, const int reducedColumn,
							 const int rowIsOdd, const int columnIsOdd)
{
	double **matrix = allocateMatrix(numRow, numColumn);
    int rowCounter;
    int columnCounter;

    // init every value of the matrix to 0
    for(rowCounter = 0; rowCounter < numRow; rowCounter++)
    {
            for(columnCounter = 0; columnCounter < numColumn; columnCounter++)
            {
                   	matrix[rowCounter][columnCounter] = 0;

            }
    }

    int arrayIndex;
    int arraySize = reducedRow * reducedColumn;
    int row, col;

    // rebuild original matrix using symmetry
    for(arrayIndex = 0; arrayIndex < arraySize; arrayIndex++)
    {
    	row = (arrayIndex/reducedColumn) + 1;
    	col = (arrayIndex % reducedColumn) + 1;

    	matrix[row][col] = array[arrayIndex];
    	matrix[row][numColumn - col - 1] = array[arrayIndex];
    	matrix[numRow - row -1][col] = array[arrayIndex];
    	matrix[numRow - row -1][numColumn - col - 1] = array[arrayIndex];
    }
	return matrix;
}

void freeArray(double *matrix)
{
	free(matrix);
}

// initialize array representing the reduced matrix with value of initial temperatures
void initArray(double *array, const int numRow, const int numColumn, const int reducedRow, const int reducedColumn)
{
	const int arraySize = ( (reducedRow * reducedColumn));
	int index;
	int rowCounter = 1;
	int columnCounter = 1;

	for(index = 0; index < arraySize; index++)
	{
		array[index] = ((columnCounter * (numColumn - columnCounter - 1)) *
				        (rowCounter * (numRow - rowCounter - 1)));

		columnCounter++;

		if(columnCounter > reducedColumn)
		{
			columnCounter = 1;
			rowCounter ++;
		}
	}
}


void calculateEulerFormulaParallel(double *array, const int reducedRow, const int reducedColumn, const int numTimeStep,
						   const double discredTime, const double subdivisionSize,
						   const int rowIsOdd, const int columnIsOdd,
						   int *aggregationArray, int rank, int numtasks)
{
	double *tempArray = NULL;
	double *previousArray = NULL;
	double *newArray = NULL;

	const double tdOverPowH = (discredTime / pow(subdivisionSize, 2));
	const double tdOverPowHTimesMinusTree = (1- (4 * tdOverPowH));
	int rowCounter;
	int columnCounter;
	int arrayIndex;

	int indexStartValue = aggregationArray[rank];
	int indexEndValue;

	if (rank < numtasks-1)
		indexEndValue = aggregationArray[rank+1] - 1;
	else
		indexEndValue = reducedRow * reducedColumn - 1;


	int rowCounterStartValue = indexStartValue / reducedColumn + 1;
	int columnCounterStartValue = indexStartValue % reducedColumn + 1;

	double above, below, left, right;
	int timeStepCounter;

	int elementPerProc = aggregationArray[1]-aggregationArray[0];	// last proc may treat less element
	const int arraySize = elementPerProc * numtasks;
	const int nbElement = reducedRow * reducedColumn;

	// allocate memory for temporary matrix
	previousArray = allocateArray(elementPerProc, numtasks);
	newArray = allocateArray(elementPerProc, numtasks);

	// copy data into temp array
	memcpy(newArray, array, sizeof(double) * arraySize);

	// for every timestep
	for(timeStepCounter = 0; timeStepCounter < numTimeStep; timeStepCounter++)
	{
		// those keep track of the index of the full matrix...not reduced
		rowCounter = rowCounterStartValue;
		columnCounter = columnCounterStartValue;

		// we swap array pointers (newArray become the previousArray)
		tempArray = previousArray;
		previousArray = newArray;
		newArray = tempArray;

		// for every element treated by this processor
		for(arrayIndex = indexStartValue; arrayIndex <= indexEndValue; arrayIndex++)
		{
			// we get the top neighbor
			if(rowCounter == 1)				{above = 0;}
			else 							{above = previousArray[arrayIndex-reducedColumn];}

			// we get the bottom neighbor
			if(rowCounter == reducedRow)
			{
				// we manage the axis
				if(rowIsOdd)						{below = previousArray[arrayIndex-reducedColumn];}
				else								{below = previousArray[arrayIndex];}
			}
			else 								{below = previousArray[arrayIndex+reducedColumn];}


			// we get the left neighbor
			if(columnCounter == 1)				{left = 0;}
			else 								{left = previousArray[arrayIndex-1];}

			// we get the right neighbor
			if(columnCounter == reducedColumn)
			{
				// we manage the axis
				if(columnIsOdd)						{right = previousArray[arrayIndex-1];}
				else								{right = previousArray[arrayIndex];}
			}
			else 								{right = previousArray[arrayIndex+1];}

			//we compute the new value
			usleep(WAIT_TIME);
			newArray[arrayIndex] = 	(tdOverPowHTimesMinusTree * previousArray[arrayIndex]) + (tdOverPowH * (above + below + left + right));

			// we increment the counters that keep track of position relatively to the full matrix
			columnCounter++;
			if(columnCounter > reducedColumn)
			{
				columnCounter = 1;
				rowCounter ++;
			}
		}
		// update with neighbor
		updateNeighbor(newArray, aggregationArray, rank, numtasks, reducedColumn);
	}
	// gather data from all the processors
	MPI_Gather (&newArray[aggregationArray[rank]], elementPerProc, MPI_DOUBLE, newArray, elementPerProc , MPI_DOUBLE, SOURCE, MPI_COMM_WORLD);

	memcpy(array, newArray, sizeof(double) * nbElement);

	freeArray(previousArray);
	freeArray(newArray);
}

void calculateEulerFormula(double *array, const int reducedRow, const int reducedColumn, const int numTimeStep,
						   const double discredTime, const double subdivisionSize,
						   const int rowIsOdd, const int columnIsOdd)
{
	double *tempArray = NULL;
	double *previousArray = NULL;
	double *newArray = NULL;

	const double tdOverPowH = (discredTime / pow(subdivisionSize, 2));
	const double tdOverPowHTimesMinusTree = (1- (4 * tdOverPowH));
	int rowCounter;
	int columnCounter;
	int arrayIndex;
	double above, below, left, right;

	const int arraySize = reducedRow * reducedColumn;
	int timeStepCounter;

	previousArray = allocateArray(reducedRow, reducedColumn);
	newArray = allocateArray(reducedRow, reducedColumn);

	memcpy(newArray, array, sizeof(double) * arraySize);


	for(timeStepCounter = 0; timeStepCounter < numTimeStep; timeStepCounter++)
	{
		// those keep track of the index of the full matrix...not reduced
		rowCounter = 1;
		columnCounter = 1;

		// we swap array pointers (newArray become the previousArray)
		tempArray = previousArray;
		previousArray = newArray;
		newArray = tempArray;

		for(arrayIndex = 0; arrayIndex < arraySize; arrayIndex++)
		{
			// we get the top neighbor
			if(rowCounter == 1)				{above = 0;}
			else 							{above = previousArray[arrayIndex-reducedColumn];}

			// we get the bottom neighbor
			if(rowCounter == reducedRow)
			{
				// we manage the axis
				if(rowIsOdd)						{below = previousArray[arrayIndex-reducedColumn];}
				else								{below = previousArray[arrayIndex];}
			}
			else 								{below = previousArray[arrayIndex+reducedColumn];}


			// we get the left neighbor
			if(columnCounter == 1)				{left = 0;}
			else 								{left = previousArray[arrayIndex-1];}

			// we get the right neighbor
			if(columnCounter == reducedColumn)
			{
				// we manage the axis
				if(columnIsOdd)						{right = previousArray[arrayIndex-1];}
				else								{right = previousArray[arrayIndex];}
			}
			else 								{right = previousArray[arrayIndex+1];}


			//we compute the new value
			usleep(WAIT_TIME);
			newArray[arrayIndex] = 	(tdOverPowHTimesMinusTree * previousArray[arrayIndex]) + (tdOverPowH * (above + below + left + right));

			columnCounter++;
			if(columnCounter > reducedColumn)
			{
				columnCounter = 1;
				rowCounter ++;
			}
		}
	}

	freeArray(previousArray);

	memcpy(array, newArray, sizeof(double) * arraySize);
}

// print the matrix properly according to weather its been rotated or not at the beginning
void printMatrix(double **matrix, const int numRows, const int numColumns, const int isReversed)
{
        int rowCounter;
        int columnCounter;

        if(isReversed)
        {
            for(columnCounter = numColumns-1; columnCounter >= 0; columnCounter--)
            {
                    for(rowCounter = 0; rowCounter < numRows; rowCounter++)
                    {
                            //The number represent the number of digits to show
                            printf(" %#4.1f", matrix[rowCounter][columnCounter]);
                    }
                    printf("\n");
            }
        }
        else
        {
            for( rowCounter = 0; rowCounter < numRows; rowCounter++)
            {
                    for(columnCounter = 0; columnCounter < numColumns; columnCounter++)
                    {
                            //The number represent the number of digits to show
                            printf(" %#4.1f", matrix[rowCounter][columnCounter]);
                    }
                    printf("\n");
            }
        }
        printf("\n");
}
